* Set working directory

cd ""

* Monte carlo analyses
do HanmerKalkan_Rep

* Create a matrix to store results
matrix results = J(3, 65, .)
matrix colnames results = cor marg_obsx2_true marg_meanx2_true marg_obsx3_true ///
marg_meanx3_true marg_obsx2_negx1 marg_meanx2_negx1 marg_obsx3_negx1 marg_meanx3_negx1 ///
mab_obsx2_negx1 mab_obsx3_negx1 mab_meanx2_negx1 mab_meanx3_negx1 rmse_obsx2_negx1 ///
rmse_obsx3_negx1 rmse_meanx2_negx1 rmse_meanx3_negx1 cov_obsx2_negx1 cov_obsx3_negx1 ///
cov_meanx2_negx1 cov_meanx3_negx1 power_obsx2_negx1 power_obsx3_negx1 power_meanx2_negx1 ///
power_meanx3_negx1 marg_obsx3_negx2 marg_meanx3_negx2 mab_obsx3_negx2 mab_meanx3_negx2 ///
rmse_obsx3_negx2 rmse_meanx3_negx2 cov_obsx3_negx2 cov_meanx3_negx2 power_obsx3_negx2 ///
power_meanx3_negx2 marg_obsx2_negx3 marg_meanx2_negx3 mab_obsx2_negx3 mab_meanx2_negx3 ///
rmse_obsx2_negx3 rmse_meanx2_negx3 cov_obsx2_negx3 cov_meanx2_negx3 power_obsx2_negx3 ///
power_meanx2_negx3 sd_obsx2_true sd_meanx2_true sd_obsx3_true sd_meanx3_true ///
sd_obsx2_negx1 sd_meanx2_negx1 sd_obsx3_negx1 sd_meanx3_negx1 sd_obsx3_negx2 ///
sd_meanx3_negx2 sd_obsx2_negx3 sd_meanx2_negx3 ///
se_obsx2_negx1 se_obsx3_negx1 se_meanx2_negx1 se_meanx3_negx1 se_obsx3_negx2 ///
se_meanx3_negx2 se_obsx2_negx3 se_meanx2_negx3

* Create a counter to store results in the matrix
local i = 1

foreach rho in 0 0.5 0.8 {
	simulate marg_obsx2_true = r(marg_obsx2_true) marg_meanx2_true = r(marg_meanx2_true) ///
	marg_obsx3_true = r(marg_obsx3_true) marg_meanx3_true = r(marg_meanx3_true) ///
	marg_obsx2_negx1 = r(marg_obsx2_negx1) marg_meanx2_negx1 = r(marg_meanx2_negx1) ///
	marg_obsx3_negx1 = r(marg_obsx3_negx1) marg_meanx3_negx1 = r(marg_meanx3_negx1) ///
	mab_obsx2_negx1 = r(mab_obsx2_negx1) mab_obsx3_negx1 = r(mab_obsx3_negx1) ///
	mab_meanx2_negx1 = r(mab_meanx2_negx1) mab_meanx3_negx1 = r(mab_meanx3_negx1) ///
	rmse_obsx2_negx1 = r(rmse_obsx2_negx1) rmse_obsx3_negx1 = r(rmse_obsx3_negx1) ///
	rmse_meanx2_negx1 = r(rmse_meanx2_negx1) rmse_meanx3_negx1 = r(rmse_meanx3_negx1) ///
	cov_obsx2_negx1 = r(cov_obsx2_negx1) cov_obsx3_negx1 = r(cov_obsx3_negx1) ///
	cov_meanx2_negx1 = r(cov_meanx2_negx1) cov_meanx3_negx1 = r(cov_meanx3_negx1) ///
	power_obsx2_negx1 = r(power_obsx2_negx1) power_obsx3_negx1 = r(power_obsx3_negx1) ///
	power_meanx2_negx1 = r(power_meanx2_negx1) power_meanx3_negx1 = r(power_meanx3_negx1) ///
	marg_obsx3_negx2 = r(marg_obsx3_negx2) marg_meanx3_negx2 = r(marg_meanx3_negx2) ///
	mab_obsx3_negx2 = r(mab_obsx3_negx2) mab_meanx3_negx2 = r(mab_meanx3_negx2) ///
	rmse_obsx3_negx2 = r(rmse_obsx3_negx2) rmse_meanx3_negx2 = r(rmse_meanx3_negx2) ///
	cov_obsx3_negx2 = r(cov_obsx3_negx2) cov_meanx3_negx2 = r(cov_meanx3_negx2) ///
	power_obsx3_negx2 = r(power_obsx3_negx2) power_meanx3_negx2 = r(power_meanx3_negx2) ///
	marg_obsx2_negx3 = r(marg_obsx2_negx3) marg_meanx2_negx3 = r(marg_meanx2_negx3) ///
	mab_obsx2_negx3 = r(mab_obsx2_negx3) mab_meanx2_negx3 = r(mab_meanx2_negx3) ///
	rmse_obsx2_negx3 = r(rmse_obsx2_negx3) rmse_meanx2_negx3 = r(rmse_meanx2_negx3) ///
	cov_obsx2_negx3 = r(cov_obsx2_negx3) cov_meanx2_negx3 = r(cov_meanx2_negx3) ///
	power_obsx2_negx3 = r(power_obsx2_negx3) power_meanx2_negx3 = r(power_meanx2_negx3) ///
	se_obsx2_negx1 = r(se_obsx2_negx1) se_obsx3_negx1 = r(se_obsx3_negx1) se_meanx2_negx1 = r(se_meanx2_negx1) ///
	se_meanx3_negx1 = r(se_meanx3_negx1) se_obsx3_negx2 = r(se_obsx3_negx2) se_meanx3_negx2 = r(se_meanx3_negx2) ///
	se_obsx2_negx3 = r(se_obsx2_negx3) se_meanx2_negx3 = r(se_meanx2_negx3), ///
	seed(99) reps(1000): hanmerkalkan `rho'
	
	matrix results[`i', 1] = `rho'
	
	* True model marginal effects
	qui su marg_obsx2_true
	matrix results[`i', 2] =r(mean)
	matrix results[`i', 46] =r(sd)
	qui su marg_meanx2_true
	matrix results[`i', 3] =r(mean)
	matrix results[`i', 47] =r(sd)
	qui su marg_obsx3_true
	matrix results[`i', 4] =r(mean)
	matrix results[`i', 48] =r(sd)
	qui su marg_meanx3_true
	matrix results[`i', 5] =r(mean)
	matrix results[`i', 49] =r(sd)
	
	
	
	* Excluded x1 model marginal effects
	qui su marg_obsx2_negx1
	matrix results[`i', 6] =r(mean)
	matrix results[`i', 50] =r(sd)
	qui su marg_meanx2_negx1
	matrix results[`i', 7] =r(mean)
	matrix results[`i', 51] =r(sd)
	qui su marg_obsx3_negx1
	matrix results[`i', 8] =r(mean)
	matrix results[`i', 52] =r(sd)
	qui su marg_meanx3_negx1
	matrix results[`i', 9] =r(mean)
	matrix results[`i', 53] =r(sd)
	
	* Excluded x1 model MAB
	qui su mab_obsx2_negx1
	matrix results[`i', 10] =r(mean)
	qui su mab_obsx3_negx1
	matrix results[`i', 11] =r(mean)
	qui su mab_meanx2_negx1
	matrix results[`i', 12] =r(mean)
	qui su mab_meanx3_negx1
	matrix results[`i', 13] =r(mean)
	
	* Excluded x1 model RMSE
	qui su rmse_obsx2_negx1
	matrix results[`i', 14] =sqrt(r(mean))
	qui su rmse_obsx3_negx1
	matrix results[`i', 15] =sqrt(r(mean))
	qui su rmse_meanx2_negx1
	matrix results[`i', 16] =sqrt(r(mean))
	qui su rmse_meanx3_negx1
	matrix results[`i', 17] =sqrt(r(mean))
	
	* Excluded x1 model COVERAGE
	qui su cov_obsx2_negx1
	matrix results[`i', 18] =r(mean)
	qui su cov_obsx3_negx1
	matrix results[`i', 19] =r(mean)
	qui su cov_meanx2_negx1
	matrix results[`i', 20] =r(mean)
	qui su cov_meanx3_negx1
	matrix results[`i', 21] =r(mean)
	
	* Excluded x1 model POWER
	qui su power_obsx2_negx1
	matrix results[`i', 22] =r(mean)
	qui su power_obsx3_negx1
	matrix results[`i', 23] =r(mean)
	qui su power_meanx2_negx1
	matrix results[`i', 24] =r(mean)
	qui su power_meanx3_negx1
	matrix results[`i', 25] =r(mean)
	
	* Excluded x2 model marginal effects
	qui su marg_obsx3_negx2
	matrix results[`i', 26] =r(mean)
	matrix results[`i', 54] =r(sd)
	qui su marg_meanx3_negx2
	matrix results[`i', 27] =r(mean)
	matrix results[`i', 55] =r(sd)
	
	* Excluded x2 model MAB
	qui su mab_obsx3_negx2
	matrix results[`i', 28] =r(mean)
	qui su mab_meanx3_negx2
	matrix results[`i', 29] =r(mean)
	
	* Excluded x2 model RMSE
	qui su rmse_obsx3_negx2
	matrix results[`i', 30] =sqrt(r(mean))
	qui su rmse_meanx3_negx2
	matrix results[`i', 31] =sqrt(r(mean))
	
	* Excluded x2 model COVERAGE
	qui su cov_obsx3_negx2
	matrix results[`i', 32] =r(mean)
	qui su cov_meanx3_negx2
	matrix results[`i', 33] =r(mean)
	
	* Excluded x2 model POWER
	qui su power_obsx3_negx2
	matrix results[`i', 34] =r(mean)
	qui su power_meanx3_negx2
	matrix results[`i', 35] =r(mean)
	
	* Excluded x3 model marginal effects
	qui su marg_obsx2_negx3
	matrix results[`i', 36] =r(mean)
	matrix results[`i', 56] =r(sd)
	qui su marg_meanx2_negx3
	matrix results[`i', 37] =r(mean)
	matrix results[`i', 57] =r(sd)
	
	* Excluded x3 model MAB
	qui su mab_obsx2_negx3
	matrix results[`i', 38] =r(mean)
	qui su mab_meanx2_negx3
	matrix results[`i', 39] =r(mean)
	
	* Excluded x1 model RMSE
	qui su rmse_obsx2_negx3
	matrix results[`i', 40] =sqrt(r(mean))
	qui su rmse_meanx2_negx3
	matrix results[`i', 41] =sqrt(r(mean))
	
	* Excluded x1 model COVERAGE
	qui su cov_obsx2_negx3
	matrix results[`i', 42] =r(mean)
	qui su cov_meanx2_negx3
	matrix results[`i', 43] =r(mean)
	
	* Excluded x1 model POWER
	qui su power_obsx2_negx3
	matrix results[`i', 44] =r(mean)
	qui su power_meanx2_negx3
	matrix results[`i', 45] =r(mean)
	
	* Mean standard errors
	qui su se_obsx2_negx1
	matrix results[`i', 58] =r(mean)
	qui su se_obsx3_negx1
	matrix results[`i', 59] =r(mean)
	qui su se_meanx2_negx1
	matrix results[`i', 60] =r(mean)
	qui su se_meanx3_negx1
	matrix results[`i', 61] =r(mean)
	qui su se_obsx3_negx2
	matrix results[`i', 62] =r(mean)
	qui su se_meanx3_negx2
	matrix results[`i', 63] =r(mean)
	qui su se_obsx2_negx3
	matrix results[`i', 64] =r(mean)
	qui su se_meanx2_negx3
	matrix results[`i', 65] =r(mean)
	
	* Iterate the counter
	local i = `i' + 1
}

* Store MC matrix results as a dataframe
clear
set obs 3
svmat results, names(col)
save hanmerkalkan_results.dta, replace

* Bias
gen bias_obsx2_negx1 = marg_obsx2_negx1 - marg_obsx2_true
gen bias_meanx2_negx1 = marg_meanx2_negx1 - marg_meanx2_true
gen bias_obsx3_negx1 =  marg_obsx3_negx1 - marg_obsx3_true
gen bias_meanx3_negx1 =  marg_meanx3_negx1 -  marg_meanx3_true
gen bias_obsx3_negx2 =  marg_obsx3_negx2 - marg_obsx3_true
gen bias_meanx3_negx2 =  marg_meanx3_negx2 - marg_meanx3_true
gen bias_obsx2_negx3 =   marg_obsx2_negx3 - marg_obsx2_true
gen bias_meanx2_negx3 =  marg_meanx2_negx3 - marg_meanx2_true

* Optimism

gen opt_obsx2_negx1 = sd_obsx2_negx1 / se_obsx2_negx1 
gen opt_obsx3_negx1 = sd_obsx3_negx1 / se_obsx3_negx1 
gen opt_meanx2_negx1 = sd_meanx2_negx1 / se_meanx2_negx1
gen opt_meanx3_negx1 = sd_meanx3_negx1 / se_meanx3_negx1 
gen opt_obsx3_negx2 = sd_obsx3_negx2 / se_obsx3_negx2 
gen opt_meanx3_negx2 = sd_meanx3_negx2 / se_meanx3_negx2
gen opt_obsx2_negx3 = sd_obsx2_negx3 / se_obsx2_negx3
gen opt_meanx2_negx3 = sd_meanx2_negx3 / se_meanx2_negx3

* Save data
save hanmerkalkan_results.dta, replace

* Reshape data
drop mab* se*
reshape long marg_ bias_ rmse_ cov_ power_ sd_ opt_, i(cor) j(effect) string
split effect, p(_) generate(newvar)
drop effect
rename newvar2 model
gen variable = substr(newvar1, -2, 2)
rename newvar1 effect
replace effect = "AME" if effect == "obsx2" | effect == "obsx3"
replace effect = "MEM" if effect == "meanx2" | effect == "meanx3"
replace model = "True Values" if model == "true"
replace model = "Model 1" if model == "negx1"
replace model = "Model 2" if model == "negx2"
replace model = "Model 3" if model == "negx3"
gen model_desc = "Excludes x1" if model == "Model 1"
replace model_desc = "Excludes x2" if model == "Model 2"
replace model_desc = "Excludes x3" if model == "Model 3"
replace model_desc = "True Marginal Effects" if model == "True Values"
rename (marg sd_ bias_ rmse_ cov_ power_ opt_ cor model effect variable model_desc) ///
 (Marginal_Effect SD Bias RMSE Coverage Power Optimism Correlation Model Effect Variable Model_Descritpion)
order Model Model_Descritpion Effect Variable Correlation Marginal Effect Bias SD Optimism RMSE Coverage Power
sort Model Effect Variable Correlation

* Save data
save hanmerkalkan_results_tidy.dta, replace

* Model 1 data
preserve
keep if Model == "Model 1" | Model == "True Values"
save HanmerKalkan_Appendix_Model1.dta, replace
restore

* Model 2 data
preserve
keep if Model == "Model 2" | Model == "True Values"
save HanmerKalkan_Appendix_Model2.dta, replace
restore

* Model 3 data
preserve
keep if Model == "Model 3" | Model == "True Values"
save HanmerKalkan_Appendix_Model3.dta, replace
restore
